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A moire pattern is formed when two copies of a periodic pattern are overlaid with a relative twist. We 
address the electronic structure of a twisted two-layer graphene system, showing that in its continuum 
Dirac model the moire pattern periodicity leads to moire Bloch bands. The two layers become more 
strongly coupled and the Dirac velocity crosses zero several times as the twist angle is reduced. For a 
discrete set of magic angles the velocity vanishes, the lowest moire band flattens, and the Dirac-point 
density-of-states and the counterflow conductivity are strongly enhanced. 

Low-energy electronic properties of few layer graphene (FLG) systems are known [l|-[8| to be strongly dependent on 
stacking arrangement. In bulk graphite 0° and 60° relative orientations of the individual layer honeycomb lattices 
yield rhombohedral and Bernal crystals, but other twist angles also appear in many samples [9]. Small twist angles 
are particularly abundant in epitaxial graphene layers grown on SiC[10], but exfoliated bilayers can also appear with 
a twist, and arbitrary alignments between adjacent layers can be obtained by folding a single graphene layer [ill [l2j. 

Recent advances in FLG preparation methods have attracted theoretical attention [l3l-[l8| to the intriguing elec- 
tronic properties of systems with arbitrary twist angles, usually focusing on the two-layer case. The problem is 
mathematically interesting because a bilayer forms a two-dimensional crystal only at a discrete set of commensurate 
rotation angles; for generic twist angles Bloch's theorem does not apply microscopically and direct electronic structure 
calculations are not feasible. For twist angles larger than a few degrees the two layers are electronically isolated to 
a remarkable degree, except at a small set of angles which yield low-order commensurate structures [HI [l8|. As the 
twist angles become smaller, interlayer coupling strengthens, and the quasiparticle velocity at the Dirac point begins 
to decrease 

Here we focus on the strongly coupled small twist angle regime. We construct a low energy continuum model 
Hamiltonian that consists of three terms: two single layer Dirac-Hamiltonian terms that account for the isolated 
graphene sheets, and a tunneling term that describes hopping between layers. The Dirac Hamiltonian [l9j for a layer 
rotated by an angle with respect to a fixed coordinate system is 



where v is the Dirac velocity, k is momentum magnitude measured from the layer's Dirac point, and the spinor 
Hamiltonian acts on an individual layer's sublattice degree-of-freedom. We choose the coordinate system depicted in 
Fig.l for which the decoupled bilayer Hamiltonian is \l)h(6/2)(l \ + |2)ft(— 0/2) (2| where projects onto layer i. 

We derive a continuum model for the tunneling term by assuming that the inter-layer tunneling amplitude between 
7r-orbitals is a smooth function of separation projected onto the graphene planes. A 7r-band tight-binding calculation 
then yields (see Supplementary Information) a continuum limit in which tunneling is local. We find that the amplitude 
for an electron residing on sublattice f3 in one layer to hop to sublattice a on the other layer is 
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where 






(2) 



w = £fc D /f} and <p = 27r/3. The three q/s in equation (pQ), depicted in Fig.l, are Dirac model momentum transfers that 
correspond to the three interlayer hopping processes. For d = and a vanishing twist angle the continuum tunneling 
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matrix is 3w5 a ASpBi independent of position. By comparing with the experimentally known [20(] electronic structure 
of an AB stacked bilayer we set w « HOmeV. The spectrum is independent of d for ^ 0. (see Supplementary 
Information). In the following we therefore set d = 0. 

The continuum model hopping Hamiltonian captures the local stacking sequence of the misaligned layers. At r = 0, 
for example, only the a = A, j3 = B element of T is non-zero because we have chosen our origin at a Bernal stacking 
point. The local lattice registration changes periodically in space with a pattern of AB points at which only Ta,b 
is not zero, BA points at which only Tb,a is non-zero and AA points at which only Ta,a — Tb,b is non-zero. The 
periodicity is controlled by the reciprocal lattice vectors q& and qt r . The magnitude of the moire lattice vector is 
therefore y/3a/[2 sin(0/2)], where a is the carbon-carbon distance in graphene[2lj. Because translational symmetry in 
the continuum model is broken only by the spatially periodic hopping Hamiltonian, we can apply Bloch's theorem to 
obtain energy bands at any twist angle 0, independent of whether the underlying lattice is commensurate. We expect 
the continuum model to be accurate up to energies of ~leV and up to angles of ~ 10°. We solve numerically for the 
moire bands using the plane- wave expansion illustrated in Fig.l. Convergence is attained by truncating momentum 
space at lattice vectors of order of w/hv. The dimension of the matrix which must be diagonalized numerically is 
roughly ~ 10 0~ 2 for small (measured in degrees), compared to the ~ 10 4 0~ 2 matrix dimension of microscopic 
tight-binding models (l3l fl5j|. 

Up to a scale factor the moire bands depend on a single parameter a = w/vkg. We have evaluated the moire bands 
as a function of their Brillouin-zone momentum k for many different twist angles; results for = 5°, 1.05° and 0.5° 
are summarized in Fig. 2. For large twist angles the low energy spectrum is virtually identical to that of an isolated 
graphene sheet, except that the velocity is slightly renormalized. Large inter-layer coupling effects appear only near 
the high energy van Hove singularities discussed by Andrei and co-workers [22j. As the twist angle is reduced, the 
number of bands in a given energy window increases and the band at the Dirac point narrows. This narrowing has 
previously been expected to develop monotonically with decreasing 0. As illustrated in Fig. 2, we instead find that the 
Dirac-point velocity vanishes already at ~ 1.05°, and that the vanishing velocity is accompanied by a very fiat moire 
band which contributes a sharp peak to the Dirac-point density-of-states (DOS). At smaller twists the Dirac-point 
velocity has a non-monotonic dependence on vanishing repeatedly at the series of magic angles illustrated in Fig. 3. 

Partial insight into the origin of these behaviors can be achieved by examining the simplest limit in which the 
momentum space lattice is truncated at the first honeycomb shell. Including the sublattice degree of freedom, this 
gives rise to the Hamiltonian 

/fc k (0/2) T b T tr T t i \ 

v = Tl fc kb (-0/2) 
Tl h^-6/2) 

V Tl /* ktl (-0/2)/ 

where k is in the moire Brillouin-zone, and kj = k + q^. This Hamiltonian acts on four two-component spinors; the 
first (^o) is at a momentum near the Dirac point of one layer and the other three ^ j are at momenta near q^ and in 
the other layer. The dependence of h{6) on angle is parametrically small and can be neglected. We have numerically 
verified that this approximation reproduces the velocity with reasonable accuracy down to the first magic angle (see 
inset of Fig. 3). It is possible to demonstrate analytically (see Supplementary Information) that this Hamiltonian has 
two zero energy states at k = 0, and a Dirac velocity renormalized by 

v* _ 1 - 3a 2 
v 1 + 6a 2 

The denominator in equation (j4]) captures the contribution of the ^j's to the normalization of the wave function while 
the numerator captures their contribution to the velocity matrix elements. For small a, equation (jl]) reduces to the 
expression v*/v = 1 — 9a 2 , first obtained by Santos et al. [Hj]. The velocity vanishes at the first magic angle because 
it is in the process of changing sign. The eigenstates at the Dirac point are a coherent combination of components in 
the two layers that have velocities of opposite sign! 

The distribution of the quasiparticle velocity between the two layers implies exotic transport characteristics for 
separately contacted layers. Consider a counterflow geometry in which the currents in the two layers flow anti-parallel 
to one another. The counter-flow velocity v CF = v(l + 3a 2 )/(l + 6a 2 ) remains finite at the magic angle when the 
bands flatten and the density of states is enhanced. A Kubo formula calculation of the counterflow conductivity then 
yields (see Supplementary Information) 

= a ° (^) , (5) 

where cfq ~ e 2 e F r/7r is the conductivity of an isolated single graphene layer. As 6 is reduced from a large value towards 
1°, v* is reduced and the density-of-states is correspondingly increased. The counterflow conductivity is enhanced 



(3) 



(4) 



3 



because of an increased density of carriers, which is not accompanied by a decrease in counterflow velocity. For a 
conventional measurement in which the current in the bilayer is unidirectional v CF in expression ([5]) is replaced by v*. 
The increase in the DOS is then exactly compensated by the reduction of the renormalized velocity and the single 
layer value of the conductivity is regained. 

In summary we have formulated a continuum model description of the electronic structure of rotated graphene 
bilayer s. The resulting moire band structure can be evaluated at arbitrary twist angles, not only at commen- 
surate values. We find that the velocity at the Dirac point oscillates with twist angle, vanishing at a series of 
magic angles which give rise to large densities-of- states and large counterflow conductivities. Many properties 
of the moire bands are still not understood. For example, although we are able to explain the largest magic 
angle analytically, the pattern of magic angles at smaller values of has so far been revealed only numerically. 
Additionally the flattening of the entire lowest moire band at « 1.05° remains a puzzle. Interesting new issues 
arise when our theory is extended to graphene stacks containing three or more layers. Finally, we remark that 
electron-electron interactions neglected is this work are certain to be important at magic twist angles in neutral sys- 
tems and could give rise to counterflow superfluidity [26, 27], flat-band magnetism (28^. or other types of ordered states. 

The authors acknowledge a helfpul conversation with Gene Mele. This work was supported by Welch Foundation 
grant F1473 and by the NSF-NRI SWAN program. 
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FIG. 1: Momentum space geometry of a twisted bilayer. a) The circles represent Dirac points of the rotated graphene 
layers, separated by kg — 2if D sin(0/2) where Kr> is the magnitude of the Brillouin-zone corner wavevector for a single layer. 
Conservation of crystal momentum implies that p' = k + qb for a tunneling process in the vicinity of the plotted Dirac points, 
b) The three equivalent Dirac points in the first Brillouin zone result in three distinct hopping processes as explained in the 
Supplementary Information. Interference between hopping processes with different wavevectors captures the spatial variation 
of inter- layer coordination that defines the moire pattern. For all the three processes \qj\ = ke] however the hopping directions 
are (0, —1) for j = 1, (\/3/2, 1/2) for j = 2 and ( — VS/2, 1/2) for j = 3. Repeated hopping generates a K-space honeycomb 
lattice. The green solid line marks the moire band Wigner-Seitz cell. In a repeated zone scheme the red and black circles mark 
the Dirac points of the two layers. 
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FIG. 2: Moire Bands, a) Energy dispersion for the 14 bands closest to the Dirac point plotted along the line A — >• B — >• C —> 
D — >• A (see Fig[T|) for (left to right) = 5°, 1.05° and 0.5°. b) Density-of-states. c) Energy as a function of twist angle for the 
k = states. Band separation decreases with as also evident from (a), d) Full dispersion of the flat band at = 1.05°. 
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FIG. 3: Renormalized Dirac-point band velocity, v* vs. a 2 where a — w/vkg for 0.18° < < 1.2°. The velocity vanishes 
for ~ 1.05°, 0.5°, 0.35°, 0.24°, and 0.2°. Inset: The renormalized velocity at larger twist angles. The solid line corresponds to 
numerical results and dashed line corresponds to analytic results based on the 8-band model. 
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I. TUNNELING MATRIX 

The tunneling matrix T^, describes a process in which an electron with momentum p' residing on sublattice f3 in 
one layer hops to a momentum state k and sublattice a in the other layer. A tight binding calculation shows that 

T a(3 _ f _k±Gj_ i[Gir a -G 2 (r^-r)-Gid] r / fi N 

kp' / j q e °k+Gi,p'+G^ \P) 

G1G2 

where Q is the unit cell area, d linearly displaces one layer relative to the other, t q is the Fourier transform of the 
tunneling amplitude t(R), r is the vector connecting the two atoms within a unit cell, r A = 0, r B = r and k = K D + k 
with K D being the Dirac momentum. In equation (j6j) we have placed the origin at a point where the A sublattice of 
one layer lies above the B sublattice of the other layer when d = 0, the sums are over reciprocal lattice vectors and 
the primed vectors p' = Mp and G 7 2 = MG2 have been transformed by the rotation matrix M. Note that the crystal 
momentum is conserved by the tunneling process because t depends only on the difference between lattice positions. 
A closely related but slightly different expression appears in Ref.[18] in which we chose the origin at a honeycomb 
lattice point. The present convention is more convenient for the discussion of small rotations relative to the Bernal 
arrangement. 

The continuum model version for T is obtained by measuring wavevectors in both layers relative to their Dirac 
points and assuming that the deviations are small compared to Brillouin-zone dimensions. t(q) vanishes rapidly with q 
on a Brillouin-zone scale (this is partly because the inter-layer distance is larger than the honeycomb lattice constant). 
Therefore only values of Gi for which |K D + Gi| = |K D | contribute significantly to T. The three possible values of 
Gi in Eq.(j6]) are therefore 0, Q^ 2 \ where the latter two vectors connect a Dirac point with its equivalent first 
Brillioun zone counterparts. (See Fig.l.) Summing the three terms we find equation (1) in the main text. 

II. INDEPENDENCE OF THE SPECTRUM ON d 

Here we show that the spectrum of misaligned bilayers is independent of linear translations of one layer with respect 
to the other using a unitary transformation that makes the Hamiltonian independent of d. Consider Hq where Q is 
a momentum in the first moire Brillouin zone. With each momentum on the k-space triangular sublattice (see Fig.l) 

k = Q + nqi + mq2 

where qi = fog (1/2, \/3/2), and q2 = ko(— 1/2, y/3/2) we associate the phase 

(/>k = nG' 2 • d + mGg • d. 

The phase associated with momentum k — k^y on the other sublattice is 0k as well. In terms of the new basis states 
exp(i(/>k)|ka) the Hamiltonian Hq is d- independent. 

III. RENORMALIZED VELOCITY IN THE 8-BAND MODEL 

For twist angles > 2° the 8-band Hamiltonian Hk, given by Eq.(3) in the main text, adequately accounts for the 
low energy spectrum. The renormalized velocity v* = <9ke£|k=o follows from the spectrum e£ of the twisted bilayer. 
We find the spectrum perturbatively in k. The Hamiltonian is expressed as a sum of the k = term T~L^ and the 
k-dependent term 

/cr-k \ 
cr-k 
cr-k 
\ cr-k/ 

and solved to leading order in 



= - 
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Consider the k = term in the Hamiltonian. The wave functions of are expressed as ^ = (^o, ^i? ^2? ^3) each 
Sfrj being a two component spinor. We now assume that %^ has zero energy eigenstates and prove our assumption 



by explicitly finding these states. The zero energy eigenstates must satisfy 



Since 



the equation for the \I>o spinor is 



T 3 h~^=0 



h V = 



(7) 



(8) 



(9) 



i.e. ^0 is one of the two zero energy states ^/q 1 ^, °f the isolated layer. The two zero energy eigenstates of T-L^ 
then follow from equations (0 and (|9j). Given that I^q^I = 1 the normalization of \£ is given by 



(10) 



with I being the identity matrix. For the second term in equation ([TO]) we use the fact that h is hermitian and the 
relations 



to obtain 



V 1 



-hj/ej 



J2TjT] = 61 



(11) 



w 



and |^| 2 = 1 + 6a 2 . 

To leading order in k the energies are therefore given by the eigenvalues of 



n^ = (i\n^\j) 



1 + Qa 2 



Using familiar Pauli matrix identities and equations (|8)llj) we obtain 

w 2 ^Tjh~ 1] cr • k/i^Tt = -a 2 ^Tjcr ■ kTj = -3a 2 cr • k. 



Aside from a renormalized velocity 



< = -v 



1 - 3a 2 a)f (j) 



is therefore identical to the continuum model Hamiltonian of single layer graphene, e£ = v*k, where the renormalized 
velocity is given by equation (4) in the main text. 



IV. COUNTERFLOW CONDUCTIVITY IN THE 8-BAND MODEL 



We find the counterflow conductivity <r CF of the bilayer system. This conductivity relates a counterflow current to 
the electric fields that induce it; the latter being oppositely orientated in the two layers. We restrict the calculations 
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to twist angles > 2° for which the 8-band model is valid and to the semiclassical regime in which e F r > 1. We 
assume that only charge carriers in the linearly dispersing sector of the lowest conduction band contribute to the 
conductivity, i.e. that the Fermi momentum is much smaller than kg and that 1/ro < Hvke where To is single particle 
lifetime. 

Using the Kubo formula we find that 

a ™ = V ^ l^-l^l 2 [^{G r kfl (e F )}} 2 (12) 

k/i 

where g = 4 accounts for the spin and valley degeneracies, 



(13) 



is the x-component of the counterflow velocity operator(we set the electric fields along the x-axis), 

GIJuj) = 1 — (14) 

is the retarded Green function with \i labeling the two Dirac bands, and e£ = /iv*k is the energy dispersion of the 
twisted bilayer at small momenta. For an electron-doped system the valence band can be neglected and 

a CF * eV/(e P ) J ^r\(^\v x c ^)\ 2 (15) 
where v* is the density of states of the twisted bilayer. The vertex function 

(^fcl^cpl^fc) =v CF cos6 k , (16) 

where v CF = v(l + 3a 2 )/(l + 6a 2 ) follows directly from the previous section if we notice the sign differences between 
the counterflow velocity operator and the normal velocity operator. The counterflow conductivity given by equation 
(5) in the main text readily follows. 



